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5  INTRODUCTION 

5.1  Nature  of  the  problem 

The  potential  benefits  of  early  diagnosis  of  cancer  were  recognized  many  years  ago,  before 
soft  tissue  imaging  was  available.  This  goal  was  behind  the  first  efforts  to  apply  ultrasound 
to  the  problem  of  the  detection/diagnosis  of  breast  cancer.  Since  then,  many  investigators 
have  devoted  effort  to  this  problem  because  of  the  known  advantages  of  ultrasound:  it  is 
non-ionizing,  relatively  inexpensive  and  uses  widely  available,  portable,  equipment.  Today, 
with  the  success  of  x-ray  mammography  as  an  early  screening  tool,  there  is  still  room  for 
improved  methods,  since  there  is  disagreement  regarding  the  use  of  ionizing  radiation  for 
screening,  or  routine  exams.  Also,  the  fibrous  tissue  of  dense  breasts  gives  poor  results  in 

conventional  mammograms. 

Ultrasound  is  now  available  as  an  adjunctive  modality  in  many  breast  clinics,  where 
it  is  used  to  determine  if  masses  with  smooth  borders  are  cystic  or  solid  and  to  examine 
dense  young  breasts.  Biopsy  is  still  used  often  to  determine  malignancy.  The  cystic/solid 
determination  can  be  a  problem  in  borderline  cases,  because  of  the  tendency  of  some 
ultrasonic  systems  to  “fill  in”  the  echo  free  space  with  artifactual  echoes  and  thus  make 
a  cyst  appear  to  be  solid.  There  is  a  need  for  better  ultrasound  systems,  as  well  as  to 
have  a  more  general-purpose  imaging  modality  available  than  the  x-ray,  particularly  if  this 

modality  had  the  advantages  listed  for  ultrasound. 

Several  factors  have  been  identified  in  the  literature  as  those  contributing  to  poor  res¬ 
olution  of  ultrasound  images.  In  B-Scan  images,  finite  signal  bandwidth  of  the  ultrasonic 
transducers  is  a  major  reason  for  low  resolution  in  the  temporal  axis,  whereas  the  non 
negligible  beam  width  highly  contributes  to  that  in  the  lateral  direction  [18],  [i],  [14],  [44], 
[21].  The  resolution  also  depends  on  the  frequency  at  which  the  imaging  system  operates. 
In  addition  to  equipment  limitations,  there  are  factors  originating  from  the  nature  of  the 
tissue  being  imaged.  Phase  aberrations  and  velocity  variations  arising  from  acoustic  in- 
homogeneitv  of  tissues  are  two  of  the  important  causes,  not  only  for  low  resolution  but 
also  for  low^contrast  in  images  [39],  [11],  [6].  The  observed  ultrasonic  image  can,  therefore, 
be  considered  as  a  distorted  version  of  the  true  tissue  image,  where  the  axial  distortion  is 
dominated  by  the  pulse-echo  wavelet  of  the  imaging  system  and  the  lateral  distortion  by 

the  lateral  beam  profile. 


5.2  Background  of  previous  work 

The  problem  of  compensating  for  imaging  distortions  has  been  very  active  research  area, 
underscoring  the  importance  in  improving  diagnostic  quality  of  ultrasonic  images.  Theie 
have  been  several  approaches  to  cancel  out  the  aberration  effects,  but  there  is  no  consensus 
as  to  the  best  way  to  achieve  it  [10],  [26],  [25],  [38].  In  [10]  it  was  proposed  estimating 
the  differences  in  arrival  times  between  two  adjacent  receiving  element  locations  using 
cross  correlation  techniques;  these  results  were  used  to  modify  phasing  characteristics  of 
the  transducers  for  ensuing  scans.  In  [25],  using  an  idea  adapted  from  optics,  phasing 
characteristics  were  determined  using  speckle  brightness  as  a  measure  of  image  quality. 

The  availability  of  convolution  models,  such  as  in  [17],  [9],  for  ultrasound  image  for¬ 
mation  and  the  wide  availability  of  digital  computers  has  given  an  added  importance  to 
discrete-time  deconvolution  methods,  as  a  means  of  improving  images  beyond  the  capabil¬ 
ities  of  hardware.  A  few  researchers  have  investigated  the  true  2-D  deconvolution  of  RF 
images  [9].  [18],  [32],  whereas  most  of  the  published  works  are  on  1-D  techniques  [44],  [21], 

[23],  [19],  [14]. 

Since  the  resolution  along  the  lateral  direction  was  much  worse  than  that  along  the 
temporal  axis,  a  number  of  attempts  were  focussed  on  deconvolution  of  lateral  image  lines. 
In  [44],  n  B-mode  image  was  considered  to  be  an  ensemble  of  lateral  lines  coiiesponding 
to  lateral  slices  through  the  envelope  detected  image,  at  given  times  (depth).  Observing 
that  the  point-spread  function  of  a  typical  pulse-echo  imaging  system  is  highly  elongated 
along  the  lateral  direction,  they  hypothesized  that  lateral  image  lines  can  be  approximately 
described  by  a  TD  convolution  model.  Their  model  consisted  of  two  1-D  terms:  a  signal  of 
interest  called  the  tissue  reflectance  and  a  blurring  kernel  in  the  lateral  direction  called  the 
lateral  point  spread  function.  The  latter  function  was  defined  to  be  the  laterally  varjdng 
component  of  the  2-D  point  spread  function,  whose  axial  variation  had  been  approximated 
by  a  Dirac  delta  function.  The  problem  of  resolution  enhancement  was  posed  as  one  of 
extracting  the  the  tissue  reflectance  from  the  observed  image,  assuming  a  perfect  knowledge 
of  the  lateral  point  spread  function.  Using  a  Gaussian  shaped  hypothetical  lateral  point 
spread  function,  it  was  shown  that  at  the  best  signal  to  noise  ratio  that  can  be  expected 
from  ultrasound  images,  deconvolution  will  lead  to  a  resolution  enhancement  of  no  better 
than  2.0.  The  definition  of  the  resolution  was  based  on  the  reciprocal  oi  the  effective  width 
of  the  lateral  point  spread  function.  However,  the  amount  of  improvement  was  also  reported 
to  be  dependent  on  the  exact  shape  of  the  lateral  point  spread  function.  These  figures  were 
found  to  be  in  agreement  with  the  empirical  numbers  reported  in  [14].  Several  others  have 
reported  results  on  1-D  lateral  deconvolution  [43],  [33],  [21]  .  where  in  [33]  it  was  concluded 
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that  the  computational  effort  on  lateral  deconvolution  was  wasted  because  of  the  very  low 
resolution  enhancement  they  could  obtain  at  the  expense  of  introducing  more  artifacts. 

The  drawbacks  of  the  lateral  deconvolution  techniques  discussed  above  are  the  following. 
Although  the  quantity  displayed  on  an  ultrasound  imager  is  the  envelope  of  the  received  RF 
signal,  the  image  formation  process  actually  occurs  in  the  RF  domain.  In  [32]  it  is  pointed 
out  that  in  general  the  convolutional  model  of  image  formation  in  the  RF  domain  does 
not  hold  for  envelope  detected  signals.  It  is  concluded  that  the.l-D  lateral  deconvolution 
on  envelope  succeeds  only  in  the  special  case  where  no  phase  interference  from  nearby 

reflectors  is  present. 

All  the  papers,  except  [21],  on  lateral  deconvolution  mentioned  above,  relied  on  mea¬ 
suring  the  lateral  beam  profile  to  be  used  in  computations.  To  obtain  the  true  tissue  image 
from  the  distorted  observation,  one  requires  quantitative  information  on  the  complex  beam 
shape  and  acoustic  velocity  variations  in-vivo,  which  are  impossible  to  measure  directly.  In 
highlv  simplified  situations  such  as  measurements  of  wire  targets  under  water,  it  is  possible 
to  get  that  information  reliably.  However,  in  the  case  of  in-vivo  tissue  targets,  such  infor¬ 
mation  is  generally  unavailable.  The  measurements  done  under  water  are  not  valid  with 
clinical  images,  even  when  the  imaging  system  used  is  the  same,  because  of  the  effects  of 
phase  aberrations,  nonlinearities  and  dispersive  attenuation  introduced  by  the  tissues  [21], 
[11],  [39]  ,  [6].  thus  dramatically  limiting  the  clinical  applicability  of  those  methods. 

In  [21] ,  a  line-by-line  lateral  deconvolution  technique  which  does  not  require  the  complex 
beam  shape  in  tissue  or  phase  information  on  adjacent  lateral  lines,  was  proposed.  This 
method  too,  however,  worked  on  amplitude  detected  B-scan  images.  It  hinges  on  the  key 
assumption  that  the  transfer  function  of  the  imaging  system  along  a  lateral  image  line 
can  be  approximated  by  the  smoothed  Fourier  transform  of  the  lateral  image  line  itself. 
A  convolution  relationship  between  the  envelope  of  a  lateral  point  spread  function  and  a 
slowly  varying  envelope  of  the  tissue  response  has  been  tacitly  assumed.  Thus  this  method 
is  subject  to  all  the  limitations  implied  by  the  above  assumptions. 

In  spite  of  the  fact  that  temporal  resolution  in  an  ultrasound  image  is  much  higher 
than  the  lateral  resolution,  axial  deconvolution  is  still  of  importance.  Besides  the  obvious 
advantage  in  improved  axial  resolution,  the  removal  of  the  effect  of  the  ultrasound  pulse 
echo  wavelet  (through  axial  deconvolution)  will  tend  to  make  the  appearance  of  images  more 
uniform  over  different  subjects  [16],  thus  simplifying  the  diagnosis  procedures.  As  the  shape 
of  the  pulse  echo  wavelet  changes  with  propagation  due  to  dispersive  attenuation,  a  first 
step  in  axial  deconvolution  often  involves  the  estimation  of  the  pulse  in  tissue.  Parametric 
modeling  of  speckle-only  image  lines  has  been  proposed  [37],  [20],  [16].  However,  such 
approaches  are  limited  by  the  problems  such  as  the  model  order  selection,  associated  with 
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parametric  modeling.  In  [19]  a  Kalman  filter  technique  was  applied  to  estimate  pulse  echo 
wavelets  as  well  as  to  simultaneously  improve  the  axial  resolution.  The  success  of  the 
method  was  reported  to  be  dependent  on  the  SNR  of  the  observations  and  the  accuracy  at 
which  the  observations  could  be  modeled.  A  non-parametric  approach  for  the  estimation 
of  the  pulse  was  proposed  in  [15],  where  the  minimum-phase  equivalent  of  the  pulse-echo 
wavelet  was  separated  from  the  tissue  response.  However,  quite  often  pulse-echo  wavelets 
and  lateral  kernels  are  non-minimum  phase  signals,  thus  limiting  the  generality  of  this 

approach. 

5.3  Methods  of  approach 

In  the  past  year  we  introduced  a  novel  non-parametric  framework  for  deconvolution  of 
B-scan  images  [1],  [2],  [3],  [4],  [5].  We  first  developed  a  model  for  the  rf  image,  and  then 
reconstructed  distortions  using  higher-order  statistics  of  the  measured  image  lines.  Based 
on  the  estimated  distortions  we  performed  deconvolution  of  the  corresponding  images  and 
demonstrated  that  the  resolution  of  ultrasound  images  of  tissue  mimicking  phantom  as  well 
as  human  tissue  images  was  significantly  improved.  In  the  past,  estimation  of  distortions 
has  been  carried  out  using  exclusively  second  order  statistics  (autocorrelation).  Autocor¬ 
relation,  however,  can  recover  only  the  minimum  phase  quivalent  of  the  true  distortions, 
because  it  is  blind  to  phase.  We  also  estimated  distortions  using  second  order  statistics 
of  the  image  lines.  We  showed  that  although  these  estimates  did  lead  to  resolution  im¬ 
provement,  the  amount  of  improvement  was  less  significant  that  the  one  obtained  with  the 
higher-order  statistics  based  estimates. 

Image  formation  process  in  the  RF-domain  is  described  by  a  2-D  convolutional  model, 
where  the  attenuation  of  the  pulse-echo  wavelet  and  beam  aberration  effects  can  be  indi¬ 
rectly  incorporated  [17],  [18],  [16].  Two  1-D  blurring  kernels,  corresponding  to  axial  and 
lateral  directions,  is  hypothesized  to  represent  distortions  along  respective  axes.  The  axial 
distortion  kernel  includes  the  blurring  effects  due  to  the  finite  bandwidth  of  the  transducer, 
and  dispersive  attenuation  of  the  pulse-echo  wavelet  in  tissue.  The  lateral  distortion  kernel 
represents  the  convolutional  components  of  lateral  blurring  due  to  the  complex  beam  pat¬ 
terns.  Formalizing  a  definition  for  resolution,  we  show  that  compensation  for  the  effects  of 

the  blurring  kernels  improves  the  resolution  of  the  image. 

The  proposed  method  has  the  advantage  of  being  able  to  estimate  these  kernels  at 
each  image  line,  axial  and  lateral,  thus  capturing  the  variations  within  the  image.  Since 
the  estimations  are  based  on  higher-order  statistics  [24]  of  RF-data,  the  estimated  kernels 
are  robust  to  additive  observation  noise  and  also  have  correct  phase.  To  the  best  of  our 
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knowledge,  the  method  we  proposed  [1]  is  the  first  one  to  in-vivo  estimate  the  distortion 
kernels  with  their  true  phase,  as  opposed  to  conventional  methods  that  estimate  minimum- 
phase  equivalent  of  kernels. 

6  BODY 

6.1  Modeling  the  rf  image 

During  an  ultrasonic  investigation,  a  three-dimensional  pulsed  pressure  field  is  emitted 
into  the  tissue.  The  field  interacts  with  the  tissue  and  part  of  it  is  reflected,  scattered,  and 
subsequently  received  by  the  transducer.  Under  the  assumptions  of  linear  propagation  and 
weak  scattering,  an  expression  for  the  received  pressure  field  was  derived  m  [1  (],  using  the 
first  order  Born  approximation.  Absorption  effects  were  neglected.  The  equation  has  been 
expressed  S/S  a  convolutional  model  in  the  following  form. 

y{r2,t)  =  Upe(f)  *t  .f(ri)  /v(D,r2,f)  +  (1) 

where: 

•  n,  r2  are  vectors  denoting  the  location  of  the  scatterer  and  the  transducer,  respectively; 

•  ”  and  denote  time  and  spatial  convolution,  respectively: 

•  f(ri)  originates  from  the  inhomogeneities  in  the  tissue  due  to  density  and  propagation 
velocity  perturbations  above  their  mean  levels,  giving  rise  to  the  back  scattered  signal  [tis- 
sue  response).: 

•  Vpe{t)  is  the  pulse-echo  wavelet  that  accounts  for  the  transducer  excitation  and  the  im¬ 
pulse  responses  during  emission  and  reception  of  the  pulse, 

•  /i.pe(ri,  r2,  f)is  the  modified  pulse-echo  spatial  impulse  response  that  relates  the  transducer 
geometry  and  the  spatial  extend  of  the  scattered  held.  The  computation  of  /ipe(ri,r2,f) 
is  based  on  the  approach  described  in  [41].  [35].  Convolutional  components  of  aberrations 
and  dispersive  attenuation,  which  introduce  spatially  varying  effects  to  the  process,  may 
be  incorporated  in  this  already  spatiall\  ^ aiding  kernel. 

•  io(r2,t)  represents  measurement  noise  and  the  unmodeled  dynamics  of  the  image  foima- 
tion  process. 

The  problem  of  extracting  the  tissue  response  /(ri)  from  the  observation  y{r2,t)  is  a 
deconvolution  problem.  Since  a  B-mode  image  is  a  mapping  from  the  3-D  tissue  space  to 
the  2-D  space  of  the  displa,y,  the  solution  is  not  unique  in  general.  This  non-uniqueness 
should  be  obviated  by  making  reasonable  assumptions  about  the  3-D  structure  of  the  tissues 
being  imaged  [9].  An  assumption  implicit  in  cases  where  the  deconvolution  is  done  using 
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kernels  confined  to  the  imaging  plane  is  that  all  image  features  in  the  imaging-plane  extend 
perpendicular  to  the  corresponding  plane  in  the  tissue  space,  so  as  to  make  the  height  of 
extension  the  effective  height  of  the  beam. 

The  convolution  model  of  (17)  expresses  the  fact  that  the  received  signal  at  the  trans¬ 
ducer  site  r2  is  a  result  of  linear  spatio-temporal  interaction  between  the  signal  of  interest 
/(n)  and  a  distortion  kernel.  Thus,  the  measured  signal  contains  a  distorted  version  of 
the  true  tissue  response  /(r,)-  Deconvolving  these  kernels  should  improve  image  resolution 
and  contrast.  It  should  also  remove  aberration-induced  artifacts  that  result  from  changing 
beam  profiles  inside  the  tissue. 

As  discussed  in  the  introduction,  efforts  to  carry  out  this  deconvolution  have  been 
hampered  by  the  difficulties  in  measuring  the  modified  spatial  impulse  response  of  the 
imaging  system.  Underwater  measurements  using  simplified  targets  would  not  reveal  any 
significant  changes  undergone  by  the  interrogating  beam  in  tissue. 

Our  goal  here  is  actually  to  identify  the  combination  of  Upe(f)  and  the  spatially  varying 
kernel  /ipe(ri, r2,  t),  and  subsequently  cancel  it  from  the  image  in  order  to  improve  lateral 
as  well  as  axial  resolution.  Let  us  combine  both  smoothing  kernels  v(f)  and  /i.pe(ri. r2, 
in  one  spatially  and  temporally  varying  kernel,  /dri,r2,f).  which  we  are  going  to  refer  to 
as  the  ultrasonic  system  impulse  response.  For  discrete  time,  and  for  some  fixed  transducei 
location  (17)  is  equivalent  to  the  following  two-dimensional  convolutional  model 

i  3 

where  y{h  n)  represents  the  sample  from  the  /-th  A-line  at  discrete  time  n.  The  goal  here 
is  to  identify  the  time  varying  ultrasonic  system  response  h{Ln),  and  recover  the  tissue 

response,  from  the  noisy  measurement  y(kj). 

However,  we  will  not  attempt  a  true  2-D  deconvolution  in  this  paper.  As  is  commonly 
done  in  ultrasound  deconvolution  literature  [15],  [23],  [18],  [19],  we  assume  that  an  RF 
A-line  can  be  expressed  as  a  convolution  between  two  1-D  axial  terms:  a  hypothetical 
tissue  response  and  a  distortion  kernel.  This  view  is  not  unique  to  the  image  deconvolution 
literature;  ultrasound  Doppler  systems  and  tissue  attenuation  estimation  techniques  tacitly 
depend  on  it  [15].  Reducing  the  problem  to  a  1-D  deconvolution  is  analogous  to  the  original 
decomposition  of  the  true  3-D  deconvolution  problem  in  to  a  2-D  one. 
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6.1.1  In  vivo  distortion  estimation  using  higher-order  statistics  (HOS) 

In  the  following  we  will  treat  each  image  line  (either  in  ateral  or  axial  direction)  separately 
assuming  the  1-D  model 

yi{n)  =  hi{n)  *  fi{n)  iUi{n),  i  =  1,2,...,  (3) 

where  z  is  the  A-line  index;  /,(n)  is  the  axial  or  lateral  tissue  response;  /z,(n)  is  the  axial  or 
lateral  distortion  kernel;  which  describes  the  distortion  associated  with  the  ?-th  line;  and  n 
denotes  discrete  time.  We  assumed  that: 


(Al)  hi{n)  is  deterministic,  possibly  non-minimum  phase, 

(A2)  fi{n)  is  stationary,  white,  independent  identically  distributed  (i.i.d.),  zero-mean  non- 
Gaussian, 

(A3)  Wi{n)  is  white  zero- mean  Gaussian,  and  independent  of 

If  we  consider  a  region  of  the  image  that  contains  only  speckle,  these  assumptions  should 
hold  reasonably  well. 

The  third  order  cumulant  of  the  observation  yi{n)  is  defined  as  [24]  . 

Cy,{T,p)  =  E{yi{n)yi(n  +  T)yi{n  +  p)]  (4) 

The  bispectrum  of  the  yi{n)  is  defined  to  be  the  Fourier  transform  of  the  third  oidei 
cumulant.  Under  assumptions  (Al)  -  (A3),  the  bispectrum  of  yi{n)  is  given  by  [24]: 


Cy,iLJl,L02)  =  UJ2) 


(5) 


where  Cf^{uJi,U2)  is  the  bispectrum  of  /,(n),  Hi{u:)  is  the  spectrum  of  hi{n)  and  (u-’i , 0)2 ) 
is  the  noise  bispectrum.  If  the  additive  noise  is  zero- mean  Gaussian,  then  Cu, ,(u.'i,a;2)  =  0 
[24]  and  (4)  becomes, 

Cy,  ( u-'i ,  tD2  )  =  C7,  ( .  ^’2 )  H,  M  H,  {uJ2)H-{^'i+u;2)  ( 6 ) 

The  bicepstrum  is  defined  as  the  cepstrum  of  the  bispectrum,  from  which  we  get: 

(777.1, m2)  =  ^7/.(777i.  777  2  )  +  C/,.  (mi ,  777  2)  , 

where  are  the  bicepstra  of  /,(.,)  and  respectively. 

If  f  in)  is  stationary  independent  identically  distributed  (i.i.d.),  its  thud  or  del  spectiuni 
is  fiat  and  equal  to  the  skewness,  v,.  of  the  process.  Therefore,  its  bicepstrum  will  be  an 
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impulse  located  at  the  origin.  In  that  case,  using  bicepstral  values  along  the  mam  axes 
except  at  the  origin,  we  can  reconstruct  a  scaled  and  shifted  version  of  hi{n)  as: 


where 


f  by,{m,0) 

I  0 


777.  >  0 

777  =  0 


(9) 


[  63^.  (-777,  0)  777  <  0. 

We  assume  that  a  similar  1-D  model  holds  in  the  lateral  direction  of  the  RF-image. 
Then,  a  similar  procedure  can  be  followed  to  estimate  lateral  distortion  kernels  at  each  line. 
The  kernels  thus  estimated  will  include  the  convolutional  components  of  aberration  as  well. 
Therefore,  axial  deconvolution  followed  by  lateral  deconvolution,  or  various  combinations 
thereof,  should  give  us  distortion  compensated  RF  images  which  will  have  higher  resolution 
and  contrast.  The  ability  to  estimate  and  remove  beam  distortion  effects  (due  to  aberration) 
is  seen  as  a  major  advantage  of  this  method  over  other  non-parametric  techniques  [15],  [21]. 
The  method  proposed  in  [15],  as  is,  can  not  estimate  non-minimum  phase  signals,  while 


that  of  [21]  has  been  designed  for  envelope  detected  signals. 


6.1.2  In  vivo  estimation  of  the  minimum  phase  equivalent  of  distortion  using 
second-order  statistics  (SOS) 

Using  the  model  of  eq.  (3),  second  order  statistics  were  also  used  to  estimate  imaging 
distortions  [12],  [13]. 

Transforming  (3)  in  the  autocorrelation  domain  (second-order  statistic)  we  get. 

r^_(r)  =  E{yiin)y,(n  +  r)}  =  72’  E  bi{n)h,{n  +  r)  +  72^‘<^(^),  (10) 

n 

where  is  the  variance  of  fi{n).  is  the  variance  of  the  noise,  and  6(r)  is  the  unit 
impulse.  The  power  cepstrum  is  defined  as  the  inverse  Fourier  transform  of  the  logaiithm 
of  the  Fourier  transform  of  the  autocorrelation.  Assuming  that  the  noise  level  is  low  enough, 
in  the  power  cepstrum  domain  we  get: 

y,{k)  ^  h{k),  ky^O.,  (11) 


where  hik)  is  the  power  cepstrum  of  the  distortion  kernel  associated  with  the  ?  th  line. 
From  its  power  cepstrum,  we  can  reconstruct  the  minimum  phase  equivalent  of  Ein)  [27] 

as :  / 1  o  \ 

/,,.(„)  =  F-Me.rp{F[ydAdu.’(A’)])],  (12) 
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where  F[.]  and  F~^[.]  denote  forward  and  inverse  Fourier  transforms,  and 


w{k)  = 


1,  k>0 

0,  otherwise. 


(13) 


6.2  Experiments 

The  goals  of  our  experimental  work  were  to: 

1.  Compare  the  axial  distortions,  estimated  based  on  higher-order  statistics  (HOS-based 
distortions),  with  experimentally  obtained  ones.  Validate  the  imposeed  assumptions 
assumptions  (see  Section  6.1.1). 

2.  Estimate  HOS-based  distortions  from  clinical  ultrasound  images.  Deconvolve  ul¬ 
trasound  images  using  the  HOS-based  estimated  distortions.  Quantify  resolution 

improvement. 

3.  Deconvolve  clinical  image  using  distortions  estimated  based  on  second-order  statistics 
(SOS-based  distortions)  as  outlined  in  Section  6.1.2).  Quantify  resolution  impiove- 

ment. 

4.  Compare  resolution  improvement  of  deconvolution  with  HOS-based  and  SOS-based 
distortions. 

Towards  the  above  goals  we  collected  the  following  data: 

•  Water-tank  measurements: 

Experiment  A 

We  obtained  B-Scan  images  of  an  ATS  model  532  contrast  resolution  phantom 
which  was  positioned  in  a  water  tank.  The  target  area  consisted  ol  the  tissue 
mimicking  background  which  had  a  scatterer  density  of  32  scatterers/irmi  V  The 
transducer  we  used  was  a  model  GE1046S2  curved  device  with  a  nominal  center 
frequency  of  3.5  MHz.  The  nominal  focal  zone  of  the  transducer  was  6-13  cm.  A 
stepper  motor  with  a  step  size  0.0257nn?  controlled  the  position  of  the  transducer 
under  the  guidance  of  a  personal  computer.  Data  acquisition  was  done  through 
a  LeCroy  model  945()A  dual  channel  digital  oscilloscope  connected  to  the  PC  by 
a  GPIB  interface.  Transducer  was  moved  across  the  scanning  plane  in  steps  of 
0.25177771.  and  RE  echos  were  sampled  at  a  13.3  MHz  rate. 
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Experiment  B 

A  long  piece  of  freshly  peeled  wire  of  diameter  0.812m7?7.  was  placed  underwater 
so  that  it  was  parallel  to  the  transducer  surface  and  perpendicular  to  the  scan¬ 
ning  plane.  It  was  kept  inside  the  focal  zone  of  the  transducer,  at  a  distance 
8cm  away  from  the  transducer  surface.  The  experimental  setup  used  here  is 
the  same  as  in  experiment  (A).  The  wavelet  reflected  from  the  wire  surface  was 
recorded  and  taken  to  be  an  approximate  estimate  of  the  pulse-echo  wavelet  of 
the  imaging  system. 

•  Measurements  from  Clinical  Equipment; 

Experiment  C 

To  demonstrate  the  performance  of  our  method  on  data  from  more  realistic 
equipment,  we  imaged  the  same  ATS532  phantom  using  a  linear  array  sector 
scan  transducer  on  a  model  UltraMark-9  clinical  imaging  system  manufactured 
by  Advanced  Technology  Laboratories,  Seattle,  U.S.A.  The  target  area  consisted 
of  two  cylinders  with  scatterer  densities  4  and  8  scatterers/mm^,  embedded 
in  a  tissue  mimicking  background  of  32  scatterers/mmA  The  nominal  center 
frequency  of  the  scanner  was  3.5  MHz:  the  field  of  vision  was  60°.  Data  were 
sampled  at  a  rate  of  12  MHz.  No  TGC  was  employed. 

Experiment  D 

To  demonstrate  the  performance  of  our  method  under  clinical  conditions,  we 
obtained  a  liver  scan  of  a  patient  imaged  by  ultrasonologists  at  the  Thomas  Jel- 
ferson  Universit}’  Hospital,  Philadelphia,  on  the  same  imaging  system  described 
under  experiment  (C).  The  patient  has  been  diagnosed  with  hypoechoic  multiple 
liver  metastatic  tumor.  No  T.G.C.  had  been  applied;  focus  of  the  imaging  sys¬ 
tem  at  transmission  was  set  at  2c?n,  and  dynamic  focusing  was  employed  with 
the  receiving  mode. 

6.3  Results 

6.3.1  Water-tank  experiments  and  the  accuracy  of  our  estimation  method 

In  order  to  demonstrate  the  validity  of  our  kernel  estimation  procedure,  we  used  the  data 
from  underwater  measurements,  i.e,  experimeni^  (A)  and  (B).  A  part  of  the  RF  image  data 
that  we  gathered  using  the  circular  transducer  in  cxpenmeni  (A)  is  shown  in  Fig.  1(a). 
where  the  logarithm  of  the  envelope  has  been  used  lor  display  purposes. 
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Under  the  assumption  that  the  axial  blurring  kernel  of  an  image  line  does  not  signifi¬ 
cantly  depend  on  the  lateral  location  of  the  line,  data  from  several  nearby  axial  RF  lines 
can  be  used  to  make  a  longer  data,  vector,  which  will  enable  us  to  obtain  better  cumulant 
estimates.  In  order  to  minimize  the  effect  of  attenuation  on  our  estimations,  from  each 
A-line  i  we  considered  data  segments  yi{k)  of  length  not  more  than  2N  samples.  The 
number  of  adjacent  line  segments  yi{k)  concatenated  was  in  the  range  i.  =  1  -  10.  In  the 

axial  kernel  estimations,  we  used  M  =  10  and  N  =  64. 

Axial  blurring  kernels  estimated  with  the  method  of  Section  6.1.1  from  different  regions 
of  the  image  in  Fig  1(a)  are  shown  in  Fig.  2(a),  in  dotted  lines.  All  of  these  kernels  have 
been  estimated  from  axial  data  obtained  at  the  same  depths  of  the  image,  but  at  different 
lateral  locations.  The  mean  axial  blurring  kernel,  ),  which  was  computed  as  the  average 
of  estimated  kernels,  is  indicated  by  the  solid  line  while  the  measured  pulse  echo  wavelet 
is  indicated  b}'  the  dashed  line. 

It  can  be  seen  that  all  of  the  estimated  kernels  possess  a  similar  structure,  resembling 
a  typical  ultrasound  pulse-echo  wavelet.  This  is  to  be  expected  since  the  component  due 
to  the  pulse-echo-wavelet  (tVe(O)  dominates  the  axial  blurring  kernel  [18].  The  variation 
among  the  estimated  kernels  may  be  attributed  to  the  statistical  estimation  errors,  de¬ 
viation  of  the  scatterer  response  from  a  statistically  white  response  (which  violates  our 
assumption  (A2))  and,  contributions  from  the  effects  of  the  medium  such  as  aberration 
and  dispersive  attenuation  which  have  spatially  varying  characteristics. 

There  is  a  reasonable  match  between  the  measured  and  average  estimated  kernels  (Fig. 
2).  within  the  limits  of  the  accuracy  of  the  experimentation.  Looking  at  their  spectra  (see 
Fig.  2(b))  we  can  clearly  see  that  the  frequency  spectrum  of  the  estimated  mean  kernel 
has  the  salient  features  contained  in  the  experimental  pulse  echo  wavelet  measured  in 
experiment  (B).  The  center  frequency  of  the  spectra  fall  very  close  to  the  nominal  center 
frequency  of  the  transducer,  3.5  MHz.  The  measured  wavelet  has  a  slightly  narrower 
main  lobe,  probably  due  to  the  ringing  introduced  by  the  wire  target  used  in  experiment 
(B).  However,  it  should  be  kept  in  mind  that  the  measured  kernel  was  obtained  from 
simplified  underwater  experiments,  which  do  not  truly  represent  the  situation  inside  the 

tissue  mimicking  phantom. 

The  differences  between  the  average  estimated  kernel  and  the  measured  kernel  can  be 
attributed  to; 

(El)  the  variance  associated  with  the  estimation  of  third  order  cumulants. 

(E2)  the  much  higher  dispersive  attenuation  and  aberrations  encountered  inside  the  phan¬ 
tom  in  experiment  (A),  compared  with  those  in  experiment  (B). 
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(E3)  the  non-whiteness  of  the  underlying  scatter  response  of  the  phantom  that  violates 
assumption  (A2), 

(E4)  approximate  realization  of  a  point  target  by  a  line  target,  in  experiment  (B). 

Averaging  techniques  that  we  used  with  both  the  cuniulant  estimates  of  observations  and 
estimated  kernels  themselves,  make  the  contributions  to  total  error  from  (El)  small.  The 
errors  introduced  by  approximating  a  point  target  by  a  thin  wire  target  does  not  introduce 
serious  errors  either;  this  method  is  being  used  in  the  testing  of  ultrasound  transducers  in 
research  environments,  [32],  [21]. 

The  major  contribution  to  the  differences  in  estimated  and  measured  kernels  are  from 
(E2)  and  (£3).  In  fact,  by  virtue  of  modifying  the  frequency  spectrum  of  the  received 
signal,  attenuations  appear  as  one  component  contributing  to  the  non-whiteness  (color)  of 
the  scatterer  response,  and  thereby  affects  the  kernel  estimation  procedure.  To  minimize 
the  effect  of  attenuations  on  the  stationarity  of  our  observations,  we  considered  only  short 
segments  (eg.  64  or  128  samples  at  a  13.3MHz  sampling  rate)  from  each  A-hne  of  inter¬ 
est.  Identifying  and  compensating  for  the  color  of  the  tissue  response  can  be  achieved  by 
applying  the  blind  deconvolution  procedure  proposed  in  [31]  on  different  image  lines  [2]. 

The  results  we  obtained  support  our  model  assumptions  (Al)-(  A3).  It  should  be  noted 
that  (A2)  actually  contradicts  the  Gaussianity  assumption,  which  is  commonly  made  for  the 
tissue  response.  However,  non  Gaussian  models  for  the  backscattered  ultrasonic  signals  have 
been  suggested  and  used  in  the  past  [42],  [34],  [29],  [30].  In  the  method  proposed  here,  if  the 
tissue  response  were  Gaussian,  it  would  have  been  suppressed  in  the  bispectrum  domain, 
thereby  rendering  the  estimation  of  an  axial  blurring  kernel  an  impossibility.  Hence,  our 
success  in  the  estimation  supports  the  non  Gaussianity  assumption  for  RF  ultrasound 
data.  The  structure  of  the  measured  kernel  {experiment  (B))  suggest  that  the  axial  kernel 
is  indeed  a  non-minimum  phase  signal;  the  method  proposed  in  this  paper  is  capable  of 
estimating  non-minimum  phase  signals  in  contrast  to  existing  techniques  which  can  only 
estimate  the  minimum-phase  equivalent  of  the  true  non-minimum  phase  signal. 

We  consider  the  agreement  between  the  estimated  and  measured  kernels  satisfactory 
within  the  abtliiy  of  underwater  measurements  (Measurement  B)  to  match  the  situation 

inside  the  phantom. 

6.3.2  HOS-based  distortion  estimation  and  the  deconvolution  of  B-mode  im¬ 
ages 

Images  of  a  tissue  mimicking  phantom  obtained  with  a  single  element  transducer. 

Having  established  the  validity  of  our  assumptions  and  the  kernel  estimation  technique 
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in  the  axial  direction,  we  proceeded  to  estimate  lateral  distortion  kernels  from  the  image 
obtained  in  experiment  (A),  and  to  perform  deconvolution.  The  lateral  distortion  kernels 
we  estimated  at  locations  covering  a  range  of  depths  6.58cm-9.42cm  from  the  transducer 
surface  can  be  seen  in  Fig.  3(a);  their  spectra  are  shown  in  Fig.  3(b).  In  these  estimations, 
we  used  N  =  32  and  M  =  10.  Five  adjacent  lateral  lines  were  used  in  each  kernel  estimation. 

In  the  lateral  direction,  the  estimated  kernels  have  characteristics  similar  to  the  lateral 
beam  profiles  of  the  imaging  system.  Unfortunately,  there  is  no  direct  way  to  verify  the 
results  of  our  lateral  estimations.  In  conventional  transverse  beam  profile  estimations,  the 
experiments  are  usually  done  under  water,  and  the  pressure  profiles  are  peak  detected.  This 
process  masks  instantaneous  features  of  the  beam,  hence  the  results  of  such  experiments 
can  not  be  used  in  a  verification.  Moreover,  it  has  been  shown  that  m  the  presence  of 
aberration,  the  lateral  point  spread  function  may  undergo  significant  modifications  [6], 
[39],  possibly  serious  enough  to  change  diagnoses  in  clinical  situations  [39].  Time  histories 
of  two  dimensional  instantaneous  beam  profiles  at  the  focus,  shown  in  [22],  aie  in  agreement 
with  the  observation  that  the  lateral  beam  profile  can  depart  heavily  from  the  ideal  in  the 
presence  of  tissue  inhomogeneities.  These  results  suggest  that,  in  principle,  underw'ater 
experiments  or  theoretical  formulae  that  do  not  take  aberration  in  to  account,  can  not  be 
used  in  verifying  lateral  blurring  kernels  estimated  Irom  complex  targets.  Further,  it  should 
be  noted  that  even  though  our  lateral  kernels  are  estimated  at  individual  lateial  image  lines 
(i.e.,  at  a  fixed  time)  nearby  lines  contribute  to  the  data  at  a  given  time  by  virtue  of  the 
spatio-temporal  nature  of  Ultrasonic  S3'stem  Impulse  Response  defined  earlier. 

Once  the  distortion  kernels  have  been  identified,  retrieving  the  corresponding  true  tis¬ 
sue  response  becomes  a  typical  deconvolution  process.  For  the  deconvolution,  we  used  the 
constrained  Wiener  filter  technique  described  in  [40].  In  performing  the  axial  deconvolu¬ 
tion,  w^e  used  the  mean  HOS-based  axial  distortion  kernel  u„(t)  shown  in  Fig.  2.,  wdiich 
was  obtained  with  the  method  of  Section  6.1.1;  in  the  lateral  direction,  the  mean  lateral 
distortion  kernel  obtained  as  the  average  of  the  estimated  kernels  shown  m  Fig.  3  was  used. 

Fig.  1(b)  shows  the  logarithmically  compressed  envelope  of  the  B-scan  image,  derived 
from  the  laterally  deconvolved  RF  data  corresponding  to  the  original  image  shown  in  Fig. 
1(a).  The  axially  deconvolved  image  is  in  Fig.  1(c)  and  the  laterally  iollow'ed  by  axially 
deconvolved  image  is  in  Fig.  1(d).  Each  image  represents  an  area  of  2.5c777X2.0cm  in  the 
true  tissue  space.  For  the  sake  of  computational  simplicity,  we  assumed  that  the  lateral 
(and  axial)  blurring  kernel  do  not  vary  significantly  over  different  RF-image  lines  or  different 
image  depths  in  the  amounts  w^e  are  concerned  with  here.  This  allowed  us  to  use  a  single 
blurring  kernel,  in  each  of  the  axial  and  lateral  directions,  to  deconvolve  the  entire  image. 
Using  the  average  lateral  kernel  to  laterall}-  deconvolve  a  region  of  the  image  is  justifiable 
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by  the  fact  that  all  the  kernels  estimated  from  the  depth  range  covering  the  image  shown 
here  (6.58cm  to  9.42c77r)  show  a  certain  degree  of  similarity.  This  may  be  due  to  the  fact 
that  the  transducer  we  used  had  a  focal  zone  of  6-13c7?7,  into  which  range  the  lateral  image 
lines  considered  above  belong.  Moreover,  the  phantom  we  used  in  experiments  had  not 
been  designed  to  simulate  strong  aberrating  effects.  .4s  for  the  axial  deconvolutions,  the 
pulse  echo  wavelet  (the  major  contributor  to  the  axial  distortion  kernel)  of  the  imaging 
system  is  known  to  stay  fairly  constant  over  different  RF-lines  across  the  image  [18],  [16]; 
our  results  shown  in  Fig.  2  also  support  that  hypothesis.  However,  the  method  proposed 
here  is  still  valid,  if  one  needs  to  estimate  kernels  for  localized  regions  or  individual  image 
lines.  To  capture  fine  details  in  the  pass  band  of  the  spectrum,  one  may  still  require  a  line 
by  line,  or  region  by  region,  lateral  deconvolution. 

Quantifying  resolution  imyrovemeni 

According  to  Fig.  1,  the  deconvolution  in  the  RF-domain  results  in  a  significant  reduc¬ 
tion  in  the  size  of  speckles,  suggesting  a  gain  in  resolution.  In  order  to  quantify  the  apparent 
increase  in  resolution,  we  defined  a  measure  of  resolution  based  on  the  2-D  anto-covariance 
function  of  the  image.  Auto-covariance  function  of  an  image  was  computed  on  the  RF  data 
corresponding  to  the  image,  with  the  peak  value  of  autocovariance  normalized  to  1.0.  The 
lateral  slice,  Loo,  through  the  peak  of  the  2-D  auto-covariance  function  was  considered  in 
defining  the  lateral  resolution.  The  axial  slice.  Aoo-  through  the  peak  in  defining  the  axral 
resolution.  Similarly,  the  lateral  resolution,  was  defined  as  the  reciprocal  of  the  width 
of  Loo.  at  d  dB  below  the  peak  of  the  slice.  Axial  resolution  was  defined  similarly  except 
for  the  fact  that  the  envelope  of  the  absolute  value  of  the  slice  Aqo  was  used.  In  this  paper 
we  used  d  -  5.00  and  d  =  10.00.  Defining  the  resolution  at  two  different  dB  levels  allows 
us  to  get  a  better  idea  of  the  shape  of  the  auto-covariance  function. 

The  reason  for  defining  the  axial  resolution  based  on  the  envelope  of  absolute  value  of 
Aoo  is  as  follows.  We  define  the  resolution  in  terms  of  the  width,  L,  of  the  mam  lobe  of  the 
auto-covariance  function  because.  L  is  a  measure  of  the  average  "smallness’  of  the  basic 
building  elements  of  the  image.  In  the  axial  direction,  the  re]u-esentation  of  a  point  target  is 
an  oscillatory,  time-limited  signal.  Therefore,  even  when  the  target  is  a  point  in  space,  the 
received  signal  would  have  an  auto-covariance  function  which  shows  oscillatory  behavior 
(eg.  Fig.  4(d)).  In  this  case,  the  central  lobe  conveying  information  on  the  average  image 
element  size  is  the  envelope  of  the  absolute  value  of  the  auto-covariance  function. 

The  resolution  was  defined  in  the  RF  domain,  because  our  kernel  estimation  and  decon¬ 
volution  procedure  was  done  entirely  in  the  RF-domain.  The  process  of  envelope  detecting 
itself  introduces  blurring  leading  to  an  inaccurate  estimate  of  the  potential  resolution  deliv- 
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ered  by  deconvolution,  thus  making  the  envelope  detected  signal  domain  a  less  than  ideal 
place  to  define  resolution.  Fig.  4(a)  shows  a  2D  shaded  mesh  plot  of  the  auto-covariance 
function  of  the  RF-data  corresponding  to  original  image  in  Fig.  1(a).  Fig.  4(b)  shows 
the  auto- variance  function  after  lateral  and  axial  deconvolutions.  Absolute  values  of  the 
auto-covariance  function  have  been  plotted  for  easy  visualization.  The  slice  Loo  is  shown  in 
Fig.  4(c),  where  the  solid  line  indicates  the  slice  corresponding  to  the  original  image  (Fig. 
1(a))  and  the  dotted  line  that  of  the  laterally  and  axially  deconvolved  image  (Fig.  1(c)). 

Corresponding  plots  of  the  axial  slice  Aqo  ^re  in  Fig.  4(d). 

The  general  decrease,  in  all  directions,  of  the  width  of  the  main-lobe  of  the  auto- 
covariance  function  due  to  deconvolution  is  obvious  from  Figs.  4(a)  and  4(b).  From  Fig. 
4(c),  we  computed  the  latcvcil  rtsolution  gain,  defined  as. 

d  =  5dB,10dB,  (14) 

where  and  respectively  represent  lateral  resolutions  before  and  after  deconvolu¬ 

tion.  A  similar  definition  holds  for  the  axial  resolution  gain.  The  lateral  resolution  gam  at 
d  =  5dB  and  d  =  lOdB  levels,  between  the  original  (Fig.  4(a))  and  the  laterally  followed 
by  axially  deconvolved  image  (Fig.  4(c))  was  found  to  be  2.7  and  3.0  respectively.  The 
corresponding  figures  for  axial  resolution  gain,  was  1.73  and  1.72  respectively. 

Our  results  indicate  that  the  lateral  resolution  gain  is  higher  than  the  axial  resolu¬ 
tion  gain.  Coupled  with  the  fact  that,  in  the  original  image  blurring  was  much  higher  in 
the  lateral  direction,  we  can  conclude  lateral  deconvolution  is  mostly  responsible  for  the 
improvement  in  overall  resolution. 

To  investigate  the  effects  of  deconvolution  on  “speckle  noise”  levels,  we  defined  the 
signal-to-noise  ratio  as,  SNR  =  (^),  where  g  is  the  mean  of  the  image  and  is  the  variance. 
Since  we  performed  our  deconvolutions  in  the  RF  domain,  our  SB R  calculations  were  also 
done  in  the  RF-domain,  on  absolute  values  of  RF-data.  Computed  over  the  original  image 
(Fig.  1(a)),  SNR  =  1.24:  computed  over  the  laterally  and  axially  deconvolved  image, 
SNR  =  1.26.  .A.fter  only  an  axial  deconvolution.  (Fig.  1(c)),  SNR  =  1.25,  and  after  only  a 
lateral  deconvolution  (Fig.  1(b)),  SNR  =  1.23.  Based  on  these  numbers,  we  conclude  that 
the  deconvolution  results  in  a  gain  in  resolution,  but  does  not  significantly  alter  speckle 
noise  levels.  The  results  of  SNR  computations  and  resolution  gains  have  been  summarized 

in  Tables  1  and  2. 

Images  of  a  tissue  mimicking  phantom  obtained  with  clinical  equipment_ 

In  order  to  investigate  the  performance  of  our  technique  with  B-scan  images  taken  from 
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modern  clinical  equipment,  we  estimated  axial  and  lateral  distortion  kernels  from  the  RF- 
data  collected  in  expenment  (C).  Fig.  5(a)  shows  the  logarithmically  compressed  envelope 
of  the  original  B-scan  image.  Average  axial  and  lateral  kernels  estimated  from  the  RF-data 
corresponding  to  Fig.  5(a)  are  illustrated  in  Fig.  6(a)  and  6(b)  respectively.  These  kernels 
indicate  the  average  of  20  kernels  estimated  from  the  tissue  mimicking  background  region 
of  the  phantom,  between  the  two  target  cylinders  (see  Fig.  5).  The  spectra  of  the  average 
axial  kernel  is  shown  in  Fig.  6(c)  and  that  of  the  average  lateral  kernel  in  Fig.  6(d). 

The  result  of  lateral  deconvolution  is  shown  in  Fig.  5(b);  the  result  of  lateral  followed 
by  axial  deconvolution  is  shown  in  Fig.  5(c).  Average  estimated  kernels  shown  in  Fig.  6 
were  used  in  both  axial  and  lateral  deconvolutions.  Clearly,  the  deconvolution  has  resulted 
in  a  reduced  speckle  size  and  cleaner,  better  defined  boundaries  of  the  target  cylinders.  It 
is  also  evident  that  the  attenuations  associated  with  the  imaging  process  show  up  much 
clearly  in  the  deconvolved  image.  This  may  have  significance  in  clinical  imaging  situations, 
where  attenuation  properties  of  tissue  convey  important  diagnostic  information. 

Fig.  7(a)  and  7(b)  illustrate  the  absolute  value  of  the  auto-covariance  function  of  the 
RF-drta  corresponding  to  the  original  (Fig.  7(a))  image  and  the  laterally  and  axially 
deconvolved  image  (Fig.  7(d)),  respectively.  Fig.  7(c)  and  7(d)  show  the  lateral  and 
axial  slices  Too  and  Aqo  used  to  compute  resolution  gains.  Based  on  data  corresponding  to 
Figs.  5  and  7,  speckle  SNR  and  resolution  gains  were  computed  and  tabulated  in  Tables 
1  and  2.  Figs.  5,  7  and  Tables  1  and  2  lead  to  the  conclusion  that  considerable  resolution 
enhancement  is  possible  with  deconvolution,  and  that  the  process  of  deconvolution  does  not 
affect  the  speckle  noise  significantly.  Once  again,  lateral  resolution  is  found  to  be  mostly 
responsible  for  the  overall  improvement  in  the  image. 

Clinical  images  of  a  human  liver 

To  evaluate  our  method  with  clinical  images,  we  obtained  a  B-scan  of  a  liver  image  as 
described  in  experiment  (D).  section  3.1.  Figure  8(a)  shows  a  part  of  the  logarithmically 
compressed  envelope  of  the  original  image.  Fig  9(a)  and  9(b)  show  the  average  axial  and 
lateral  distortions  estimated  from  the  RF-data  corresponding  to  Fig.  8(a).  Their  spectia 
can  be  seen  in  Figs.  9(c)  and  9(d). 

The  results  of  lateral  deconvolution  is  shown  in  Fig.  8(b):  results  of  lateral  followed  by 
axial  is  in  Fig.  8(c).  In  all  cases,  logarithmically  compressed  envelope  of  the  images  has 

been  displayed. 

Again,  to  visualize  the  improvement  in  resolution,  we  plotted  the  auto-covariances  of 
the  original  and  deconvolved  images.  Fig.  10(a)  and  10(b)  illustrates  the  shaded  2-D  mesh 
plot  of  the  auto-covariances  of  RF-data  corresponding  to  Figs.  8(a)  and  8(c)  respectively. 
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Corresponding  Loo  and  Aoo  slices  are  shown  in  Figs.  10(c)  and  10(d).  Slices  corresponding 
to  original  image  are  shown  in  solid  lines,  whereas  those  of  the  deconvolved  image  aie  shown 

in  dotted  lines. 

The  lateral  resolution  gain,  at  d  =  bdB  and  d  =  lOdB  was  found  to  be  3.0  and  2.5 
respectively  Corresponding  numbers  for  the  axial  resolution  gain  were  1.7  and  1.5.  Speckle 
SNR  ratio^omputed  for  the  original  image  (Fig.  8(a))  was  SNR  =  1.20;  after  lateral  and 
axial  deconvolutions  (Fig.  8(c))  SNR  =  1.33.  Tables  1  and  2  summarize  these  results. 
These  results  suggest  that  considerable  gain  in  resolution  is  possible  with  axial  and 

lateral  deconvolution  of  clinical  images. 

6.3.3  SOS-based  estimated  distortions  and  deconvolution  of  B-mode  images 

From  the  images  obtained  at  experiments  A,  C  and  D.  SOS-based  distortions  were  estimated 
with  the  method  of  Section  6.1.2.  The  images  were  subsequently  deconvolved,  and  the 
deconvolution  results  are  shown  in  Figs.  11,  12  and  13,  respectively.  The  corresponding 
resolution  gains  and  speckle  SNR  are  shown  in  Tables  3  and  4. 

Comparing  Tables  1  and  2  to  3  and  4,  we  can  conclude  that  the  resolution  improvement, 
both  axial  and  lateral,  was  consistently  superior  when  deconvolution  was  performed  with 
the  HOS-based  estimated  distortions.  Deconvolution  with  HOS-based  estimates  led  to 
axial  resolution  improvement  1.15  -  1.9  times  higher  than  deconvolution  with  SOS-based 
estimates.  The  improvement  of  lateral  resolution  when  deconvolving  using  HOS  based 
estimates  was  more  dramatic;  it  was  better  than  the  SOS  based  deconvolution  by  a  factor 
in  the  range  1.68-4.5.  The  superiority  of  the  HOS-based  deconvolution  is  also  evident 
by  comparing  the  corresponding  images  deconvolved  with  HOS  and  SOS-based  distortion 

estimates. 

T  CONCLUSIONS/FUTURE  WORK 

Processing  ultrasound  images  with  higher-order  spectral  operations  we  were  able  to  identify 
the  distortion  introduced  by  the  ultrasonic  system  and  the  medium.  The  proposed  method 
makes  it  possible  to  estimate  both  axial  and  instantaneous  lateral  blurring  kernels,  working 
on  B-mode  RF  data.  Distortion  identification  and  sulisequent  cancelling  (decon\olution) 
operations  were  carried  out  on  1-D  lines  of  the  RF  image,  thereby  obviating  the  theoretical 
difficulties  faced  liy  earlier  attempts  at  beam  deconvolution  on  envelope  detected  images. 
The  method  is  capable  of  estimating  mixed  phase  distortion  kernels,  and  is  immune  to 

additive  Gaussian  noise. 
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Performing  underwater  experiments  employing  single  element  transducers,  commercial 
tissue  mimicking  phantoms  and  simulated  point  targets,  we  showed  that  our  kernel  estima¬ 
tion  procedure  could  be  done  with  reasonable  accuracy.  The  accuracy  of  the  estimations 

were  verified  by  measuring  axial  kernels  undei\\atei. 

Deconvolution  results  obtained  with  phantom  data  and  clinical  images  indicate  con¬ 
siderable  resolution  improvement.  Lateral  deconvolution  contributes  heavily  to  the  gam 
in  resolution,  where  the  resolution  was  defined  in  terms  of  the  dimensions  of  the  auto¬ 
covariance  function  of  the  image.  These  results  are  significant  because  in  unprocessed 
images,  lateral'  distortions  are  known  to  be  as  much  as  3-5  times  severe  than  axial  distor¬ 
tions,  leading  to  a  change  in  appearance  in  clinical  images  depending  on  the  orientation  of 
the  scanner.  Through  lateral  deconvolution,  one  can  compensate  for  the  lateral  distortions 
and  try  to  achieve  consistent  images  independent  from  the  angular  position  of  the  scanner. 
In  the  past,  distortion  estimation  was  carried  out  exclusively  in  the  second-order  statistics 
domain.  We  demonstrated  that  deconvolution  based  on  HOS-based  distortion  estimates 
leads  to  superior  axial  and  in  particular  lateral  resolution  improvement  than  the  one  with 
SOS-based  estimates. 

It  has  been  reported  that  high  transducer  frequencies  required  for  high  resolution  imag¬ 
ing  actually  lead  to  lower  resolution  in  the  presence  of  increased  aberration  m  breast  at 
higher  frequencies  [8].  In  [36]  it  is  pointed  out  that  medium  inhomogeneities  are  also  impor¬ 
tant  in  systems  where  a  larger  aperture  is  used  for  higher  lateral  resolution.  In  general,  the 
aberration  correction  is  of  fundamental  importance  in  high  resolution  medical  ultrasonic 
systems.  The  ability  of  the  proposed  method  to  compensate  for  convolution  components 
of  aberration  holds  promise  in  improving  the  diagnostic  value  of  B-mode  medical  images 
beyond  the  capabilities  of  hardware. 

Deconvolution,  viewed  as  an  image  de-blurring/restoration  operation,  should  reveal 
those  small  structures  (such  as  early  stage  tumors)  that  had  been  hidden  away  from 
view  due  to  imaging  distortions  (including  convolutional  components  of  aberrations),  but 
otherwise  would  show  up  on  a  hypothcHcal.  ptrjeci  imager.  A  deconvolved  image  also 
means  access  to  a  distortion-free  tissue  signal,  which  is  largely  independent  of  the  imaging 
system.  Therefore,  tissue  characterization  schemes  which  ideally  require  imagmg-system- 
independent  data  could  be  based  on  deconvolved  RF  images.  In  addition,  the  estimated 
distortion  kernels  themselves  carry  information  on  the  statistical  structure  of  scatterer  field 
as  manifested  through  the  color  of  the  scatterer  response  (tissue  response),  attenuations 
and.  propagation  non-linearities  associated  with  the  phantom  (tissue). 

Currently,  efforts  are  under  way  to  see  how  much  improvement  in  resolution  is  possible 
with  ultrasound  images  of  the  breast.  We  have  ordered  a  custom  made  anthropomorphic 
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expect  to  begin  experiments  with  breast  data  within  August 


breast  phantom  and  we  we 
1995. 
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9  FIGURES/TABLES 


Fig.l  A  speckle-only  part  of  the  ultrasound  image  of  the  tissue  mimicking  phantom  obtained 
with  a  focused  single  element  transducer  experiment  (A)-,  (b)  the  result  of  lateral  de- 
convolution;  (c)  the  result  of  axial  deconvolution  and,  (d)  the  result  of  lateral  followe 
by  axial  deconvolution.  In  all  cases,  the  logarithmically  compressed  envelope  shoj^. 
Deconvolution  was  performed  with  the  higher-order  statistics  based  estimated  dis¬ 
tortions. 


Fig. 2  Underwater  experiments  {^experiments  (A)  and  (B)):  (a)  axial  kernels  estimated  at  var¬ 
ious  lateral  positions  in  the  B-mode  image  of  the  tissue  mimicking  phantom  (dotted 
lines);  average  of  the  estimated  kernels  (solid  line).  The  experimental  kernel,  measured 
as  the  reflection  off  a  0.812mm  diameter  wire  surface,  under  water.  The  close  agreement 
between  the  estimated  and  the  measured  kernels  (with  in  limits  of  experimental  errors), 
indicates  the  success  of  the  estimation  procedure,  (b)  Spectra  of  the  average  estimated 
kernel  (solid  line)  and  the  measured  kernel  (dashed  line).  The  center  frequency  of  both 
kernels  are  compatible  with  the  nominal  center  frequency  of  the  transducer,  3.5MHz. 


Fig.3  (a)  Lateral  distortion  kernels  estimated  at  various  axial  depths  from  the  RF  data  cor 
responding  to  the  image  shown  in  Fig.  1(a).  All  the  kernels  are  estimated  within  the 
focal  zone  of  the  transducer,  6  -  13cm.  (b)  Spectra  of  the  lateral  kernels  shown  in  (a). 
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Fig.4  Shaded  auto-covariance  function  of  the  RF-data  corresponding  to:  (a)  the  original  image 
shoAvn  in  Fig.  1(a);  (b)  the  laterally  and  axially  deconvolved  image  shown  in  Fig.  1(d). 
Absolute  values  of  the  auto-covariance  functions  have  been  shown  for  easy  visualization, 
(c)  The  lateral  slice  Loo  and  (d)  the  axial  slice  Aoo  of  the  auto-covariance  function 
corresponding  to:  the  original  image  (solid  lines),  the  deconvolved  image  (dotted  lines). 


(C) 


Fig.5  {Experiment  (C)):  (a)  The  original  image  of  the  tissue  mimicking  phantom  obtained 
with  the  linear  array  transducer  on  a  clinical  imaging  system;  (b)  the  result  of  lateral 
deconvolution  and  (c)  the  result  of  lateral  and  axial  deconvolution.  The  logarithmically 
compressed  envelope  has  been  used  for  display.  Deconvolution  was  performed  with  the 
higher-order  statistics  based  estimated  distortions. 


Fig.6  {Experiment  (C)):  (a)  The  average  of  the  estimated  ax.al  kerne  s  from 

Losing  to  image  shown  in  Fig.  6(a),  and  (b)  its  spectrum,  (c)  The  ^ge  of  the 
iLral  kernels  estimated  from  RF-data  corresponding  to  Fig.  5(a)  and  (d)  its  epectrum. 
These  kernels  were  estimated  from  the  region  between  the  two  cylindrical  targets  in  F  g. 

5(a). 
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Fig.7  Shaded  auto-covariance  function  of  the  RF-data  corresponding  to:  (a)  the  origma 

shown  in  Fig.  5(a);  (b)  the  laterally  and  axially  deconvolved  image  shown  in  Fig.  5(d). 
Absolute  values  of  the  auto-covariance  functions  have  been  shown  for  easy  visualization, 
(c)  The  lateral  slice  Loo  and  (d)  the  axial  slice  Aoo  of  the  auto- covariance  function 
corresponding  to:  the  original  image  (solid  lines),  the  deconvolved  image  (dotted  lines). 


Fig.8 


(Experiment  (D)):  (a)  A  part  of  the  clinical  image  of  a  human  liver  containing  a  tumor, 
(b)  the  result  of  lateral  deconvolution  and  (c)  the  result  of  lateral  and  axial  deconvolu¬ 
tion.  In  all  cases  the  logarithm  of  the  envelope  is  shown.  Deconvolution  was  performed 
with  the  higher-order  statistics  based  estimated  distortions. 


Fig.9  (Experiment  (D)):  (a)  The  average  of  the  estimated  axial  kernels  from  RF-data  corre- 
spending  to  image  shown  in  Fig.  8(a),  and  (b)  its  spectrum,  (c)  The  average  of  the 
lateral  kernels  estimated  from  RF-data  corresponding  to  Fig.  8(a)  and  (d)  its  spectrum. 
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Fig.  10  Shaded  auto-covariance  function  of  the  RF-data  corresponding  to:  (a)  the  original 

shown  in  Fig.  8(a);  (b)  the  laterally  and  axially  deconvolved  image  shown  in  Fig.  8(d). 
Absolute  values  of  the  auto-covariance  functions  have  been  shown  for  easy  visualization, 
(c)  The  lateral  slice  Loo  and  (d)  the  axial  slice  Aoo  of  the  auto-covariance  function 
corresponding  to:  the  original  image  (solid  lines),  the  deconvolved  image  (dotted  lines). 
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Fis  11  A  speckle-only  part  of  the  ultrasound  image  of  the  tissue  mimicking  phantom  obtained 
with  a  focused  single  element  transducer  experiment  (A)\  (b)  the  result  of  lateral  de 
convolution;  (c)  the  result  of  axial  deconvolution  and,  (d)  the  result  of  lateral  followe 
by  axial  deconvolution.  In  all  cases,  the  logarithmically  comprised  envelope  ‘s  shown. 
Deconvolution  was  performed  with  the  second-order  statistics  based  estimated  dis- 

tortions. 


Fig. 12  [Experiment  (C)):  (a)  The  original  image  of  the  tissue  mimicking  phantom  obtained 
with  the  linear  array  transducer  on  a  clinical  imaging  system;  (b)  the  result  of  lateral 
deconvolution  and  (c)  the  result  of  lateral  and  axial  deconvolution.  The  logarithmically 
compressed  envelope  has  been  used  for  display.  Deconvolution  was  performed  with  the 
second-order  statistics  based  estimated  distortions. 


original  image  laterally  deconvolved  image 


laterally  and  axially  deconvolved  image 


Fig. 13  {Experiment  (D)):  (a)  A  part  of  the  clinical  image  of  a  human  liver  containing  a  tumor; 
(b)  the  result  of  lateral  deconvolution  and  (c)  the  result  of  lateral  and  axial  deconvolu¬ 
tion.  In  all  cases  the  logarithm  of  the  envelope  is  shown.  Deconvolution  was  performed 
with  the  second-order  statistics  based  estimated  distortions. 


Table  1; 


Table  2: 
ratio 


Resolution  gains  due  to  deconvolution  with  HOS-based  estimated  distrotions 


Resolution  Gains 

lateral 

axial 

5dB 

lOdB 

5dB 

lOdB 

.. 

Experiment  (A) 

2.7 

3.on 

1.7 

1.7 

Experiment  ( C) 

5.2 

4.2 

1.8 

1.9 

Experiment  (D) 

3.0 

2.5 

1.7 

1.5 

The  effect  of  deconvolution  with  HOS-based  distortions  on  the  speckle  SNR 


SNR 

original 

deconvolved 

Experiment  (A) 

1.24 

1.26 

Experiment  (C) 

1.17 

1.10 

Experiment  (D) 

1.20 

1.30 
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Table  3: 
ratio 


Table  4: 
ratio 


The  effect  of  deconvolution  with  SOS-based  distortions  on  the  speckle  SNR 


Resolution  Gains 

lateral 

axial 

5dB 

lOdB 

5dB 

lOdB 

Experiment  (A) 

1.6 

1.35 

1.53 

1.20 

Experiment  (C) 

1.85 

0.93 

1.48 

1.64 

Experiment  (D) 

1.50 

1.44 

1.50 

1.23 

The  effect  of  deconvolution  with  SOS-based  distortions  on  the  speckle  SNR 


SNR 

original 

deconvolved 

Experiment  (A) 

1.21 

1.24 

Experiment  ( C) 

1.2 

1.20 

Experiment  (D) 

1.10 

1.17 
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MEETING  ABSTRACTS 


SPIE  International  Symposium  on  Optics,  Imaging  mid  Instrumentation,  San  Diego,  July 
1994. 

Blind  Deconvolution  of  Ultrasound  Images 

Udantha  R.  Abeyratne,  Athina  P.  Petropulu*  and  John  M.  Reid 

Biomedical  Engineering  and  Science  Institute 
*Department  of  Electrical  and  Computer  Engineering, 

Drexel  University,  Philadelphia,  PA  19104 

We  address  the  problem  of  improving  the  resolution  of  ultrasound  images  using  blind 
deconvolution.  The  transducer  measurement  that  forms  the  ultrasound  image  can  be  ex¬ 
pressed  as  the  convolution  of  two  terms,  the  tissue  response  and  the  ultrasonic  system 
response,  plus  additive  noise.  The  quantity  of  interest  is  the  tissue  response,  however,  due 
to  the  convolution  operation,  we  measure  a  blurred  version  of  it,  which  obscures  the  fine 
details  in  the  image.  Our  goal  is  to  remove  the  blurring  caused  by  the  ultrasonic  system, 
in  order  to  enhance  the  diagnostic  quality  of  the  ultrasound  image.  Under  the  assump¬ 
tion  that  speckle  behaves  as  an  i.i.d.  zero-mean  non-Gaussian  process,  we  were  able  to 
reconstruct  the  blurring  kernel  using  bicepstrum  operations  on  corresponding  A-scan  lines. 
Based  on  the  estimated  blurring  kernel  we  performed  deconvolution  on  the  measured  im¬ 
age.  Preliminary  results  obtained  from  ultrasound  images  of  a  tissue  mimicking  phantom 
indicate  significant  resolution  improvement. 


20^^  International  Symposium  on  Ultrasonic  Imaging  and  Tissue  Characterization,  Ailing- 
ton,  VA,  June  1995. 

Estimating  Imaging-Distortions  and  the  Color  of  the  Tissue  Response  From 

Ultrasonic  Images 

Udantha  R.  Abeyratne,  Athina  P.  PctropuhC,  John  M.  Reid  and  Thomas  Golas* 

Biomedical  Engineering  and  Science  Institute 
“Department  of  Electrical  and  Computer  Engineering, 

Drexel  University,  Philadelphia,  PA  19104 

We  address  the  problem  of  estimating  imaging-distortions  and  modeling  the  tissue  re¬ 
sponse  from  clinical  B-scan  ultrasound  images,  based  on  rf  A-lines  that  form  the  image.  We 


45 


model  rf  A-lines  as:  yi{n)  =  h,{n)  *  /,(n)  -f  w,{n).t  =  1,2,  •  •  •,  where  t/,(n)  is  the  observed 
rf  A-line,  fi{n)  is  the  tissue  response,  Wi{n)  is  the  observation  noise  and  hi{n)  is  the  ker¬ 
nel  representing  imaging  distortions;  the  symbol  stands  for  the  convolution  operation. 
The  fol-  lowing  assumptions  are  made:  [Al]  hi{n)  is  deterministic,  possibly  non-  minimum 
phase,  [A2]  fi{n)  is  stationary,  zero-mean  non-Gaussian,  modeled  as  the  convolution  of  a 
white  component  e,(n)  and  a  deterministic  part  t,(n)  corresponding  to  the  statistical  color 
of  the  issue  response,  [A3]  Wi{n)  is  zero-mean,  Gaussian,  and  independent  of  ft{n).  A  col 
ored  random  process  is  taken  as  the  model  for  the  tissue  response,  considering  the  fact  that 
under-  lying  scatterers  need  not  be  randomly  and  independently  located.  Applying  higher 
order  spectra  based  blind  deconvolution  on  pairs  of  closely  located  rf  A-lines,  we  identify 
the  color  of  the  tissue  response  as  well  as  the  distortion  kernels.  The  color  of  the  tissue 
response  is  independent  of  the  imaging  system  and  can  be  used  in  tissue  characterization. 
We  estimate  distortion  kernels  and  the  color  of  the  tissue  response,  from  a  B-mode  image 
of  a  human  liver,  containing  a  tumor.  The  statistical  color  of  the  tissue  response  estimated 
inside  the  tumor  has  features  that  are  distinctly  different  from  those  estimated  outside  the 
tumor.  This  may  be  attributed  to  the  change  in  the  tissue  structures  brought  about  by 
the  tumor.  Our  kernel  estimation  and  color  identification  process  is  based  on  data  fiom 
regions  as  small  as  12  mm^  (approx.)  in  the  tissue  space.  The  advantages  of  the  proposed 
approach  are:(a)  since  the  estimations  are  carried  out  in  higher  order  spectra  domain,  results 
are  immune  to  additive  Gaussian  noise  in  observations,  and  (b)  both  minimum  phase  and 
non-minimum  phase  kernels  can  be  estimated.  We  compare  our  method  with  existing  non- 
parametric  pulse  estimation  techniques  and  demonstrate  relative  strengths  and  weaknesses. 


IEEE  Wo7'kshop  on  Nonlinear  Signal  and  Image  Processing,  Neos  Marmaras,  Greece.  June 
1995. 

Higher  Order  Spectra  Based  Deconvolution  of  Ultrasound  Images 

Udantha  R.  Abejuatne,  Athina  P.  Petropulut,  m.  Reid 

Biomedical  Engineering  and  Science  Institute 
t  Electrical  and  Computer  Engineering  Department 
Drexel  University,  Philadelphia,  P.^  19104 

Our  goal  is  to  model  and  identify  the  tissue  response  based  on  the  backscattered  rf 
signal  that  forms  the  ultrasound  image.  We  model  the  rl  ultrasonic  image,  in  both  axial 
and  lateral  directions,  as  the  convolution  between  an  1-D  hypothetical  tissue  response  and 
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the  ultrasonic  system  response  (distortion),  plus  additive  Gaussian  noise.  We  model  the 
tissue  response  as  a  non-Gaussian,  colored  random  process.  Closely  spaced  axial  (lateral) 
image  lines  contain  as  a  common  convolutional  term  the  axial  (lateral)  distortion,  whereas 
the  noncommon  terms  are  due  to  the  tissue  color.  Applying  blind  deconvolution,  we  identify 
the  color  of  the  tissue  response  as  well  as  the  corresponding  distortion  kernels.  The  color 
of  the  tissue  response  is  independent  of  the  imaging  system  and  can  be  used  as  a  tissue 
characterization  parameter.  Estimated  kernels  are  compared  with  experimentally  obtained 
kernels,  and  the  deconvolution  of  real  ultrasound  images  is  performed. 
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Summary 


It  is  well  established  that  identification  of  a  nonmininium  phase  non-Gaussian  random 
process  can  be  achieved  based  on  its  higher-order  cumulants  (order  three  or  higher).  We 
propose  a  computationally  attractive  approach  for  system  identification,  whose  complex¬ 
ity  remains  constant  as  the  order  of  the  employed  cumulants  increases.  We  define  the 
Q- weighted  n'^’-order  cumulant  projection  of  the  process  x{k)  to  be 


Pn 


(T;a)  =  E-  E  ^  complex, \a\  =  1, 


:i5’) 


T2 


Tt7-1 


where  c^(-)  denotes  order  cumulant.  For  the  process 

x{k)  =  h{k)  *  e{k). 


where  e(k)  is  white  non-Gaussian  zero-mean,  and  h(k)  is  deterministic  generally  nonmini 
mum  phase,  we  prove  that 


where  P^(:‘,<y) 
shown  that  by 


P^(z:a)  =  cH(a-^=)H{o 


(71-1)  ^-1 


(17 


is  the  Z-transform  of  p^(r;a).  and  c  is  a  constant.  ;From  (17).  it  can  be 
controlling  a.  H{z)  can  be  identified  for  n  >  2.  If  p?;(/-';a)  and  h{k)  denote 
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(18) 


the  complex  cepstra  of  a)  and  k{k),  respectively,  we  show  that 


h{k) 

h{-k) 


pg(fc;a)-Q^”  A’  >  0 

Q,k  —  (yi^  ^)k 

-pl{t,a)  +  a^Pl{k-,l) 


(19) 


which  leads  to  a  scaled  and  shifted  version  of  hik)  via  inverse  cepstrum  operations.  We 
show  that  we  can  always  find  a  complex  number  a  with  unit  magnitude  to  guarantee  that 
the  denominator  in  (18)  and  (19)  is  different  from  zero,  and  discuss  the  rationale  behind 
choosing  a  particular  value  for  a. 

Based  on  equation  (17)  a  formula  for  the  estimation  of  the  Fourier  phase  of  h{k)  from 
the  phase  of  the  o-weighted  cumulant  projection  is  also  deiived. 
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